#This file is an abbreviated translation of the R code here:
    # https://r-forge.r-project.org/scm/viewvc.php/pkg/skellam/R/pskellam.sp.r?view=markup&revision=19&root=healthqueues

#Libraries
import numpy as np
import scipy.stats as stats
from mpmath import mp

#Skellam CDF
def skcdf(q, lambda1, lambda2):
    # Luganni-Rice saddlepoint CDF with Butler's 2nd continuity correction
    #Assume lower-tailed
    xm = -np.floor(q) - 0.5          #Continuity corrected x

    #Distribution-specific variables
    s = np.log(0.5*(xm + np.sqrt(xm**2 + 4*lambda2*lambda1))/lambda2) #Saddlept
    K = lambda2 * (np.exp(s) - 1) + lambda1 * (np.exp(-s) - 1)     #CGF(s)
    K2 = lambda2*np.exp(s) + lambda1*np.exp(-s)    #CGF''(s)
    
    #Code depending on distribution only through previous variables
    u2 = 2*np.sinh(0.5*s)*np.sqrt(K2)
    w2 = np.sign(s)*np.sqrt(2*(s*xm - K))
    
    #For the normal draws, we need to begin to mpmathify
    ret = mp.fsub(mp.ncdf(-w2),mp.fmul(mp.npdf(w2),(1/w2-1/u2)))
    
    #Avoid numeric problems near the removable discontinuity
    xe = (xm + (lambda1 - lambda2))/np.sqrt(lambda1 + lambda2)
    g1 = (lambda1 - lambda2)/(lambda1 + lambda2)**1.5
    ew = abs(xe) < 1e-4
    if ew:
        ret = mp.fadd(mp.ncdf(-xe),mp.fmul(mp.npdf(xe),g1/6*(1 - xe^2)))
    return ret

#Skellam PMF
def skpmf(x, lambda1, lambda2):
    #Saddlepoint density (PMF) for Skellam distribution
    s = np.log(0.5*(x + np.sqrt(x**2 + 4*lambda1*lambda2))/lambda1)#Saddlepoint
    K = lambda1*(np.exp(s) - 1) + lambda2*(np.exp(-s) - 1)	# CGF(s)
    K2 = lambda1*np.exp(s) + lambda2*np.exp(-s)        # CGF''(s)
    c = (1 - ((lambda1*np.exp(s) - lambda2*np.exp(-s))/K2)**2*5/3)/K2*0.125 + 1
    ret = np.exp(K - x*s)/np.sqrt(2*np.pi*K2)*(1 + c)*0.5
    return ret

